home *** CD-ROM | disk | FTP | other *** search
/ The Best of Down Under Games / The Best of Down Under Games.iso / 3dfx Screen Savers / Chrome / Sphere.c < prev    next >
C/C++ Source or Header  |  1997-05-20  |  9KB  |  312 lines
  1. /*% cc -g sphere.c -o sphere -lm
  2.  *
  3.  * sphere - generate a triangle mesh approximating a sphere by
  4.  *  recursive subdivision. First approximation is an platonic
  5.  *  solid; each level of refinement increases the number of
  6.  *  triangles by a factor of 4.
  7.  *
  8.  * Level 3 (128 triangles for an octahedron) is a good tradeoff if
  9.  *  gouraud shading is used to render the database.
  10.  *
  11.  * Usage: sphere [level] [-p] [-c] [-f] [-t] [-i]
  12.  *      level is an integer >= 1 setting the recursion level (default 1).
  13.  *      -c causes triangles to be generated with vertices in counterclockwise
  14.  *          order as viewed from the outside in a RHS coordinate system.
  15.  *          The default is clockwise order.
  16.  *      -t starts with a tetrahedron instead of an octahedron
  17.  *      -i starts with a icosahedron instead of an octahedron
  18.  *
  19.  *  The subroutines print_object() and print_triangle() should
  20.  *  be changed to generate whatever the desired database format is.
  21.  *
  22.  * Jon Leech (leech@cs.unc.edu) 3/24/89
  23.  * icosahedral code added by Jim Buddenhagen (jb1556@daditz.sbc.com) 5/93
  24.  */
  25. #include <stdio.h>
  26. #include <ctype.h>
  27. #include <math.h>
  28. #include <stdlib.h>
  29. #include <malloc.h>
  30.  
  31. typedef struct {
  32.     double  x, y, z;
  33. } point;
  34.  
  35. typedef struct {
  36.     point     pt[3];    /* Vertices of triangle */
  37.     double    area;     /* Unused; might be used for adaptive subdivision */
  38. } triangle;
  39.  
  40. typedef struct {
  41.     int       npoly;    /* # of triangles in object */
  42.     triangle *poly;     /* Triangles */
  43. } object;
  44.  
  45. /* Six equidistant points lying on the unit sphere */
  46. #define XPLUS {  1,  0,  0 }    /*  X */
  47. #define XMIN  { -1,  0,  0 }    /* -X */
  48. #define YPLUS {  0,  1,  0 }    /*  Y */
  49. #define YMIN  {  0, -1,  0 }    /* -Y */
  50. #define ZPLUS {  0,  0,  1 }    /*  Z */
  51. #define ZMIN  {  0,  0, -1 }    /* -Z */
  52.  
  53. /* Vertices of a unit octahedron */
  54. triangle octahedron[] = {
  55.     { { XPLUS, ZPLUS, YPLUS }, 0.0 },
  56.     { { YPLUS, ZPLUS, XMIN  }, 0.0 },
  57.     { { XMIN , ZPLUS, YMIN  }, 0.0 },
  58.     { { YMIN , ZPLUS, XPLUS }, 0.0 },
  59.     { { XPLUS, YPLUS, ZMIN  }, 0.0 },
  60.     { { YPLUS, XMIN , ZMIN  }, 0.0 },
  61.     { { XMIN , YMIN , ZMIN  }, 0.0 },
  62.     { { YMIN , XPLUS, ZMIN  }, 0.0 }
  63. };
  64.  
  65. /* A unit octahedron */
  66. object oct = {
  67.     sizeof(octahedron) / sizeof(octahedron[0]),
  68.     &octahedron[0]
  69. };
  70.  
  71. /* Vertices of a tetrahedron */
  72. #define sqrt_3  0.5773502692
  73. #define PPP {  0.5773502692,  0.5773502692,  0.5773502692 }   /* +X, +Y, +Z */
  74. #define MMP { -0.5773502692, -0.5773502692,  0.5773502692 }   /* -X, -Y, +Z */
  75. #define MPM { -0.5773502692,  0.5773502692, -0.5773502692 }   /* -X, +Y, -Z */
  76. #define PMM {  0.5773502692, -0.5773502692, -0.5773502692 }   /* +X, -Y, -Z */
  77.  
  78. /* Structure describing a tetrahedron */
  79. triangle tetrahedron[] = {
  80.     {{ PPP, MMP, MPM }, 0.0},
  81.     {{ PPP, PMM, MMP }, 0.0},
  82.     {{ MPM, MMP, PMM }, 0.0},
  83.     {{ PMM, PPP, MPM }, 0.0} };
  84.  
  85. object tet = {
  86.     sizeof(tetrahedron) / sizeof(tetrahedron[0]),
  87.     &tetrahedron[0]
  88. };
  89.  
  90. /* Twelve vertices of icosahedron on unit sphere */
  91. #define tau 0.8506508084      /* t=(1+sqrt(5))/2, tau=t/sqrt(1+t^2)  */
  92. #define one 0.5257311121      /* one=1/sqrt(1+t^2) , unit sphere     */
  93. #define ZA {  tau,  one,    0 }
  94. #define ZB { -tau,  one,    0 }
  95. #define ZC { -tau, -one,    0 }
  96. #define ZD {  tau, -one,    0 }
  97. #define YA {  one,   0 ,  tau }
  98. #define YB {  one,   0 , -tau }
  99. #define YC { -one,   0 , -tau }
  100. #define YD { -one,   0 ,  tau }
  101. #define XA {   0 ,  tau,  one }
  102. #define XB {   0 , -tau,  one }
  103. #define XC {   0 , -tau, -one }
  104. #define XD {   0 ,  tau, -one }
  105.  
  106. /* Structure for unit icosahedron */
  107. triangle icosahedron[] = {
  108.     { { YA, XA, YD }, 0.0 },
  109.     { { YA, YD, XB }, 0.0 },
  110.     { { YB, YC, XD }, 0.0 },
  111.     { { YB, XC, YC }, 0.0 },
  112.     { { ZA, YA, ZD }, 0.0 },
  113.     { { ZA, ZD, YB }, 0.0 },
  114.     { { ZC, YD, ZB }, 0.0 },
  115.     { { ZC, ZB, YC }, 0.0 },
  116.     { { XA, ZA, XD }, 0.0 },
  117.     { { XA, XD, ZB }, 0.0 },
  118.     { { XB, XC, ZD }, 0.0 },
  119.     { { XB, ZC, XC }, 0.0 },
  120.     { { XA, YA, ZA }, 0.0 },
  121.     { { XD, ZA, YB }, 0.0 },
  122.     { { YA, XB, ZD }, 0.0 },
  123.     { { YB, ZD, XC }, 0.0 },
  124.     { { YD, XA, ZB }, 0.0 },
  125.     { { YC, ZB, XD }, 0.0 },
  126.     { { YD, ZC, XB }, 0.0 },
  127.     { { YC, XC, ZC }, 0.0 }
  128. };
  129.  
  130. /* A unit icosahedron */
  131. object ico = {
  132.     sizeof(icosahedron) / sizeof(icosahedron[0]),
  133.     &icosahedron[0]
  134. };
  135.  
  136. int Flatflag = 0;   /* Don't generate per-vertex normals */
  137.  
  138. /* Forward declarations */
  139. point *normalize(/* point *p */);
  140. point *midpoint(/* point *a, point *b */);
  141. void flip_object(/* object *obj */);
  142. void print_object(/* object *obj */);
  143. void print_triangle(/* triangle *t */);
  144.  
  145. /* returns the number of triangles in the sphere */
  146. object *sphere(maxlevel)
  147. int maxlevel;
  148. {
  149.     object *old = &oct,         /* Default is octahedron */
  150.            *new;
  151.     int     ccwflag = 0,        /* Reverse vertex order if true */
  152.             i,
  153.             level;              /* Current subdivision level */
  154.  
  155.  
  156.         /* flip_object(old); */
  157.  
  158.     /* Subdivide each starting triangle (maxlevel - 1) times */
  159.     for (level = 1; level < maxlevel; level++) {
  160.         /* Allocate a new object */
  161.         new = (object *)malloc(sizeof(object));
  162.         if (new == NULL) {
  163.             fprintf(stderr, "Out of memory on subdivision level %d\n",level);
  164.             exit(1);
  165.         }
  166.         new->npoly = old->npoly * 4;
  167.  
  168.         /* Allocate 4* the number of points in the current approximation */
  169.         new->poly  = (triangle *)malloc(new->npoly * sizeof(triangle));
  170.         if (new->poly == NULL) {
  171.             fprintf(stderr, "Out of memory on subdivision level %d\n",level);
  172.             exit(1);
  173.         }
  174.  
  175.         /* Subdivide each triangle in the old approximation and normalize
  176.          *  the new points thus generated to lie on the surface of the unit
  177.          *  sphere.
  178.          * Each input triangle with vertices labelled [0,1,2] as shown
  179.          *  below will be turned into four new triangles:
  180.          *
  181.          *                      Make new points
  182.          *                          a = (0+2)/2
  183.          *                          b = (0+1)/2
  184.          *                          c = (1+2)/2
  185.          *        1
  186.          *       /\             Normalize a, b, c
  187.          *      /  \
  188.          *    b/____\ c         Construct new triangles
  189.          *    /\    /\              [0,b,a]
  190.          *   /  \  /  \             [b,1,c]
  191.          *  /____\/____\            [a,b,c]
  192.          * 0      a     2           [a,c,2]
  193.          */
  194.         for (i = 0; i < old->npoly; i++) {
  195.             triangle
  196.                  *oldt = &old->poly[i],
  197.                  *newt = &new->poly[i*4];
  198.             point a, b, c;
  199.  
  200.             a = *normalize(midpoint(&oldt->pt[0], &oldt->pt[2]));
  201.             b = *normalize(midpoint(&oldt->pt[0], &oldt->pt[1]));
  202.             c = *normalize(midpoint(&oldt->pt[1], &oldt->pt[2]));
  203.  
  204.             newt->pt[0] = oldt->pt[0];
  205.             newt->pt[1] = b;
  206.             newt->pt[2] = a;
  207.             newt++;
  208.  
  209.             newt->pt[0] = b;
  210.             newt->pt[1] = oldt->pt[1];
  211.             newt->pt[2] = c;
  212.             newt++;
  213.  
  214.             newt->pt[0] = a;
  215.             newt->pt[1] = b;
  216.             newt->pt[2] = c;
  217.             newt++;
  218.  
  219.             newt->pt[0] = a;
  220.             newt->pt[1] = c;
  221.             newt->pt[2] = oldt->pt[2];
  222.         }
  223.  
  224.         if (level > 1) {
  225.             free(old->poly);
  226.             free(old);
  227.         }
  228.  
  229.         /* Continue subdividing new triangles */
  230.         old = new;
  231.     }
  232.  
  233.     /* Print out resulting approximation */
  234.  
  235.     print_object(old);
  236.  
  237.   return(old);
  238. }
  239.  
  240. /* Normalize a point p */
  241. point *normalize(p)
  242. point *p;
  243. {
  244.     static point r;
  245.     double mag;
  246.  
  247.     r = *p;
  248.     mag = r.x * r.x + r.y * r.y + r.z * r.z;
  249.     if (mag != 0.0) {
  250.         mag = 1.0 / sqrt(mag);
  251.         r.x *= mag;
  252.         r.y *= mag;
  253.         r.z *= mag;
  254.     }
  255.     return &r;
  256. }
  257.  
  258. /* Return the midpoint on the line between two points */
  259. point *midpoint(a, b)
  260. point *a, *b;
  261. {
  262.     static point r;
  263.  
  264.     r.x = (a->x + b->x) * 0.5;
  265.     r.y = (a->y + b->y) * 0.5;
  266.     r.z = (a->z + b->z) * 0.5;
  267.     return &r;
  268. }
  269.  
  270. /* Reverse order of points in each triangle */
  271. void flip_object(obj)
  272. object *obj;
  273. {
  274.     int i;
  275.     for (i = 0; i < obj->npoly; i++) {
  276.         point tmp;
  277.                        tmp = obj->poly[i].pt[0];
  278.         obj->poly[i].pt[0] = obj->poly[i].pt[2];
  279.         obj->poly[i].pt[2] = tmp;
  280.     }
  281. }
  282.  
  283. /* Write out all triangles in an object */
  284. void print_object(obj)
  285. object *obj;
  286. {
  287.     int i;
  288.  
  289.     /* Spit out coordinates for each triangle */
  290.     for (i = 0; i < obj->npoly; i++)
  291.         print_triangle(&obj->poly[i]);
  292.  
  293. }
  294.  
  295. /* Output a triangle */
  296. void print_triangle(t)
  297. triangle *t;
  298. {
  299.  
  300.         /* Modify this to generate your favorite output format
  301.          * Triangle vertices are in t->pt[0..2].{x,y,z}
  302.          * A generic format is provided.
  303.         printf("triangle\n");
  304.         for (i = 0; i < 3; i++)
  305.             printf("\t%g %g %g\n", t->pt[i].x, t->pt[i].y, t->pt[i].z);
  306.          */
  307.     
  308. }
  309.  
  310.  
  311.  
  312.